######################################################################
# Syntax: Political Polarization during the COVID-19 Pandemic
# 
# R-Version: 3.6.1
# Author: Sebastian Jungkunz (University of Duisburg-Essen)
# E-Mail: sebastian.jungkunz@uni-due.de
######################################################################

rm(list=ls(all=TRUE))

#packages
library(tidyverse)
library(sjmisc)
library(sjPlot)
library(sjlabelled)
library(ggeffects)
library(gghighlight)
library(cowplot)
library(texreg)
library(ggpubr)
library(lubridate)
library(Hmisc)
library(dotwhisker)
library(stargazer)
library(tidyr)

select <- dplyr::select
filter <- dplyr::filter

# set working directory
setwd("SET_YOUR_WD")


## load data
d <- read_stata("./Data/corona_data_2020.dta")

#reshape data for for robustness checks
d_short <- d %>% select(id, wave, group, feel, eco_sal, vote, edu, sex, income)
d_wide <- reshape(d_short, idvar = "id", timevar = "wave", direction = "wide")



#################################################   
####### Predictions: Distress on Sympathy #######
#################################################   

### Attitudes

# Nationality

nat_ita_eco_sal_w1 <- lm(feel ~ eco_sal+ edu + sex + income, data = subset(d, (group==9 & wave==0)))
nat_ita_eco_sal_w1_p <- ggeffect(nat_ita_eco_sal_w1, terms = c("eco_sal [1:5]"), ci.lvl = 0.95)
nat_ita_eco_sal_w2 <- lm(feel ~ eco_sal+ edu + sex + income, data = subset(d, (group==9 & wave==1)))
nat_ita_eco_sal_w2_p <- ggeffect(nat_ita_eco_sal_w2, terms = c("eco_sal [1:5]"), ci.lvl = 0.95)

nat_tur_eco_sal_w1 <- lm(feel ~ eco_sal+ edu + sex + income, data = subset(d, (group==8 & wave==0)))
nat_tur_eco_sal_w1_p <- ggeffect(nat_tur_eco_sal_w1, terms = c("eco_sal [1:5]"), ci.lvl = 0.95)
nat_tur_eco_sal_w2 <- lm(feel ~ eco_sal+ edu + sex + income, data = subset(d, (group==8 & wave==1)))
nat_tur_eco_sal_w2_p <- ggeffect(nat_tur_eco_sal_w2, terms = c("eco_sal [1:5]"), ci.lvl = 0.95)


# robustness check: use wave1 IVs to predict wave2 DV
nat_ita_eco_sal_w2_w1 <- MASS::rlm(feel.1 ~ eco_sal.0 + edu.0 + sex.0 + income.0, data = subset(d_wide, (group.1==9)))
nat_ita_eco_sal_w2_w1_p <- ggeffect(nat_ita_eco_sal_w2_w1, terms = c("eco_sal.0 [1:5]"), ci.lvl = 0.95)

nat_tur_eco_sal_w2_w1 <- MASS::rlm(feel.1 ~ eco_sal.0 + edu.0 + sex.0 + income.0, data =  subset(d_wide, (group.1==8)))
nat_tur_eco_sal_w2_w1_p <- ggeffect(nat_tur_eco_sal_w2_w1, terms = c("eco_sal.0 [1:5]"), ci.lvl = 0.95)

nat_ita_eco_sal_w2_w1_p$conf.low <- nat_ita_eco_sal_w2_w1_p$predicted - 1.96 * nat_ita_eco_sal_w2_w1_p$std.error 
nat_ita_eco_sal_w2_w1_p$conf.high <- nat_ita_eco_sal_w2_w1_p$predicted + 1.96 * nat_ita_eco_sal_w2_w1_p$std.error
nat_tur_eco_sal_w2_w1_p$conf.low <- nat_tur_eco_sal_w2_w1_p$predicted - 1.96 * nat_tur_eco_sal_w2_w1_p$std.error 
nat_tur_eco_sal_w2_w1_p$conf.high <- nat_tur_eco_sal_w2_w1_p$predicted + 1.96 * nat_tur_eco_sal_w2_w1_p$std.error


# use function to calculate R2: https://stats.stackexchange.com/questions/83826/is-a-weighted-r2-in-robust-linear-model-meaningful-for-goodness-of-fit-analys
r2ww <- function(x){
  SSe <- sum(x$w*(x$resid)^2)
  observed <- x$resid+x$fitted
  SSt <- sum(x$w*(observed-weighted.mean(observed,x$w))^2)
  value <- 1-SSe/SSt;
  return(value);
}

r2ww(nat_ita_eco_sal_w2_w1)
r2ww(nat_tur_eco_sal_w2_w1)




###Plots

#nationality
nat_ita_eco_sal_w1_p$group <- 2
nat_ita_eco_sal_w2_p$group <- 2
nat_ita_eco_sal_w2_w1_p$group <- 2

nat_tur_eco_sal_w1_p$group <- 3
nat_tur_eco_sal_w2_p$group <- 3
nat_tur_eco_sal_w2_w1_p$group <- 3


nateco_sal_w1 <- rbind(nat_ita_eco_sal_w1_p, nat_tur_eco_sal_w1_p) 
nateco_sal_w1$group <- as.factor(nateco_sal_w1$group)
nateco_sal_w2 <- rbind(nat_ita_eco_sal_w2_p, nat_tur_eco_sal_w2_p) 
nateco_sal_w2$group <- as.factor(nateco_sal_w2$group)
nateco_sal_w2_w1 <- rbind(nat_ita_eco_sal_w2_w1_p, nat_tur_eco_sal_w2_w1_p) 
nateco_sal_w2_w1$group <- as.factor(nateco_sal_w2_w1$group)



#wave 1
plot_nat_eco_sal_w1 <- ggplot(nateco_sal_w1, aes(x=x, y=predicted, group=group)) +
  theme_minimal() + 
  geom_line(aes(linetype=group, color="black"), size=0.7) +
  geom_ribbon(aes(ymin=nateco_sal_w1$conf.low, ymax=nateco_sal_w1$conf.high), linetype=0, alpha=0.1, show.legend = F) +
  scale_y_continuous(limit = c(-50, 50)) +
  scale_x_continuous(limit = c(1,5), breaks=c(1,2,3,4,5)) +  
  xlab("") + 
  ylab("Sympathy") + 
  labs(color  = "Party", linetype = "Party", fill = "Party", group="Party") +
  theme(legend.position = "bottom", legend.margin=margin(0,0,0,0), legend.box.margin=margin(-5,-5,-5,-5), legend.title=element_blank(),
        legend.key.size = unit(2,"line")) + 
  scale_fill_manual(values=c("grey", "grey"), name=NULL) + 
  scale_color_manual(values=c("black", "black")) +
  scale_linetype_manual(values=c("solid", "dashed"), labels = c("Italian", "Turkish")) +
  theme(legend.position = c(0.93, 0.29),
        legend.justification=c(1, 1), 
        legend.background = element_rect(colour="grey80"),
        legend.margin=margin(-3,5,5,5),
        legend.key.size = unit(1.3, "lines")) +
  theme(legend.title=element_blank()) +
  guides(color=FALSE)
plot_nat_eco_sal_w1

ggsave("./Figures/nat_distress_w1.pdf", height = 3 , width = 5, scale=1.4, device = cairo_pdf)


#wave 2
plot_nat_eco_sal_w2 <- ggplot(nateco_sal_w2, aes(x=x, y=predicted, group=group)) +
  theme_minimal() + 
  geom_line(aes(linetype=group, color="black"), size=0.7) +
  geom_ribbon(aes(ymin=nateco_sal_w2$conf.low, ymax=nateco_sal_w2$conf.high), linetype=0, alpha=0.1, show.legend = F) +
  scale_y_continuous(limit = c(-50, 50)) +
  scale_x_continuous(limit = c(1,5), breaks=c(1,2,3,4,5)) +  
  xlab("") + 
  ylab("") + 
  labs(color  = "Party", linetype = "Party", fill = "Party", group="Party") +
  theme(legend.position = "bottom", legend.margin=margin(0,0,0,0), legend.box.margin=margin(-5,-5,-5,-5), legend.title=element_blank(),
        legend.key.size = unit(2,"line")) + 
  scale_fill_manual(values=c("grey", "grey"), name=NULL) + 
  scale_color_manual(values=c("black", "black")) +
  scale_linetype_manual(values=c("solid", "dashed"), labels = c("Italian", "Turkish")) +
  theme(legend.position = c(0.93, 0.29),
        legend.justification=c(1, 1), 
        legend.background = element_rect(colour="grey80"),
        legend.margin=margin(-3,5,5,5),
        legend.key.size = unit(1.3, "lines")) +
  theme(legend.title=element_blank()) +
  guides(color=FALSE)
plot_nat_eco_sal_w2

ggsave("./Figures/nat_distress_w2.pdf", height = 3 , width = 5, scale=1.4, device = cairo_pdf)


# robustness check
plot_nat_eco_sal_w2_w1 <- ggplot(nateco_sal_w2_w1, aes(x=x, y=predicted, group=group)) +
  theme_minimal() + 
  geom_line(aes(linetype=group, color="black"), size=0.7) +
  geom_ribbon(aes(ymin=nateco_sal_w2_w1$conf.low, ymax=nateco_sal_w2_w1$conf.high), linetype=0, alpha=0.1, show.legend = F) +
  scale_y_continuous(limit = c(-50, 50)) +
  scale_x_continuous(limit = c(1,5), breaks=c(1,2,3,4,5)) +  
  xlab("Economic Distress") + 
  ylab("Sympathy") + 
  labs(color  = "Party", linetype = "Party", fill = "Party", group="Party") +
  theme(legend.position = "bottom", legend.margin=margin(0,0,0,0), legend.box.margin=margin(-5,-5,-5,-5), legend.title=element_blank(),
        legend.key.size = unit(2,"line")) + 
  scale_fill_manual(values=c("grey", "grey"), name=NULL) + 
  scale_color_manual(values=c("black", "black")) +
  scale_linetype_manual(values=c("solid", "dashed"), labels = c("Italian", "Turkish")) +
  theme(legend.position = c(0.93, 0.29),
        legend.justification=c(1, 1), 
        legend.background = element_rect(colour="grey80"),
        legend.margin=margin(-3,5,5,5),
        legend.key.size = unit(1.3, "lines")) +
  theme(legend.title=element_blank()) +
  guides(color=FALSE)
plot_nat_eco_sal_w2_w1

ggsave("./Figures/nat_distress_w2_w1.pdf", height = 3 , width = 5, scale=1.4, device = cairo_pdf)




#### Asymmetric relationships


# afd vs. all others and vice versa
afd_main_eco_sal_w1 <- lm(feel ~ eco_sal + edu+ sex + income, data = subset(d, ((group>=2 & group<=5) & (vote==6) & wave==0)))
afd_main_eco_sal_w1_p <- ggeffect(afd_main_eco_sal_w1, terms = c("eco_sal [1:5]"), ci.lvl = 0.95)
afd_main_eco_sal_w2 <- lm(feel ~ eco_sal+ + edu+ sex + income, data = subset(d, ((group>=2 & group<=5) & (vote==6) & wave==1)))
afd_main_eco_sal_w2_p <- ggeffect(afd_main_eco_sal_w2, terms = c("eco_sal [1:5]"), ci.lvl = 0.95)

main_afd_eco_sal_w1 <- lm(feel ~ eco_sal + edu+ sex + income, data =  subset(d, (group==1 & (vote!=6 & vote!=8) & wave==0)))
main_afd_eco_sal_w1_p <- ggeffect(main_afd_eco_sal_w1, terms = c("eco_sal [1:5]"), ci.lvl = 0.95)
main_afd_eco_sal_w2 <- lm(feel ~ eco_sal + edu+ sex + income, data =  subset(d, (group==1 & (vote!=6 & vote!=8) & wave==1)))
main_afd_eco_sal_w2_p <- ggeffect(main_afd_eco_sal_w2, terms = c("eco_sal [1:5]"), ci.lvl = 0.95)




#afd vs. all others: use w1 to predict w2

afd_main_eco_sal_w2_w1 <- MASS::rlm(feel.1 ~ eco_sal.0 + edu.0 + sex.0 + income.0, data = subset(d_wide, ((group.1>=2 & group.1<=5) & (vote.1==6))))
afd_main_eco_sal_w2_w1_p <- ggeffect(afd_main_eco_sal_w2_w1, terms = c("eco_sal.0 [1:5]"), ci.lvl = 0.95)

main_afd_eco_sal_w2_w1 <- MASS::rlm(feel.1 ~ eco_sal.0 + edu.0 + sex.0 + income.0, data =  subset(d_wide, (group.1==1 & (vote.1!=6 & vote.1!=8))))
main_afd_eco_sal_w2_w1_p <- ggeffect(main_afd_eco_sal_w2_w1, terms = c("eco_sal.0 [1:5]"), ci.lvl = 0.95)

afd_main_eco_sal_w2_w1_p$conf.low <- afd_main_eco_sal_w2_w1_p$predicted - 1.96 * afd_main_eco_sal_w2_w1_p$std.error 
afd_main_eco_sal_w2_w1_p$conf.high <- afd_main_eco_sal_w2_w1_p$predicted + 1.96 * afd_main_eco_sal_w2_w1_p$std.error
main_afd_eco_sal_w2_w1_p$conf.low <- main_afd_eco_sal_w2_w1_p$predicted - 1.96 * main_afd_eco_sal_w2_w1_p$std.error 
main_afd_eco_sal_w2_w1_p$conf.high <- main_afd_eco_sal_w2_w1_p$predicted + 1.96 * main_afd_eco_sal_w2_w1_p$std.error


r2ww(afd_main_eco_sal_w2_w1)
r2ww(main_afd_eco_sal_w2_w1)



###Plots
afd_main_eco_sal_w1_p$group <- 1
afd_main_eco_sal_w2_p$group <- 1
afd_main_eco_sal_w2_w1_p$group <- 1

main_afd_eco_sal_w1_p$group <- 2
main_afd_eco_sal_w2_p$group <- 2
main_afd_eco_sal_w2_w1_p$group <- 2

#afd_lw_p$group <- 7
#lw_afd_p$group <- 8


party_afd_main_eco_sal_w1 <- rbind(afd_main_eco_sal_w1_p, main_afd_eco_sal_w1_p) 
party_afd_main_eco_sal_w1$group <- as.factor(party_afd_main_eco_sal_w1$group)
party_afd_main_eco_sal_w2 <- rbind(afd_main_eco_sal_w2_p, main_afd_eco_sal_w2_p) 
party_afd_main_eco_sal_w2$group <- as.factor(party_afd_main_eco_sal_w2$group)
party_afd_main_eco_sal_w2_w1 <- rbind(afd_main_eco_sal_w2_w1_p, main_afd_eco_sal_w2_w1_p) 
party_afd_main_eco_sal_w2_w1$group <- as.factor(party_afd_main_eco_sal_w2_w1$group)


#wave 1
plot_party_afd_main_eco_sal_w1 <- ggplot(party_afd_main_eco_sal_w1, aes(x=x, y=predicted, group=group)) +
  theme_minimal() + 
  geom_line(aes(linetype=group, color="black"), size=0.7) +
  geom_ribbon(aes(ymin=party_afd_main_eco_sal_w1$conf.low, ymax=party_afd_main_eco_sal_w1$conf.high), linetype=0, alpha=0.1, show.legend = F) +
  scale_y_continuous(limit = c(-50, 50)) +
  scale_x_continuous(limit = c(1,5), breaks=c(1,2,3,4,5)) +  
  xlab("Economic Distress") + 
  ylab("Sympathy") + 
  labs(color  = "Party", linetype = "Party", fill = "Party", group="Party") +
  theme(legend.position = "bottom", legend.margin=margin(0,0,0,0), legend.box.margin=margin(-5,-5,-5,-5), legend.title=element_blank(),
        legend.key.size = unit(2,"line")) + 
  scale_fill_manual(values=c("grey", "grey"), name=NULL) + 
  scale_color_manual(values=c("black", "black")) +
  scale_linetype_manual(values=c("solid", "dashed"), labels = c("AfD vs. Other", "Other vs. AfD")) +
  theme(legend.position = c(0.94, 0.92),
        legend.justification=c(1, 1), 
        legend.background = element_rect(colour="grey80"),
        legend.margin=margin(-3,5,5,5),
        legend.key.size = unit(1.3, "lines")) +
  theme(legend.title=element_blank()) +
  guides(color=FALSE)
plot_party_afd_main_eco_sal_w1

ggsave("./Figures/party_afd_main_distress_w1.pdf", height = 3 , width = 5, scale=1.4, device = cairo_pdf)


#wave 2
plot_party_afd_main_eco_sal_w2 <- ggplot(party_afd_main_eco_sal_w2, aes(x=x, y=predicted, group=group)) +
  theme_minimal() + 
  geom_line(aes(linetype=group, color="black"), size=0.7) +
  geom_ribbon(aes(ymin=party_afd_main_eco_sal_w2$conf.low, ymax=party_afd_main_eco_sal_w2$conf.high), linetype=0, alpha=0.1, show.legend = F) +
  scale_y_continuous(limit = c(-50, 50)) +
  scale_x_continuous(limit = c(1,5), breaks=c(1,2,3,4,5)) +  
  xlab("Economic Distress") + 
  ylab("") + 
  labs(color  = "Party", linetype = "Party", fill = "Party", group="Party") +
  theme(legend.position = "bottom", legend.margin=margin(0,0,0,0), legend.box.margin=margin(-5,-5,-5,-5), legend.title=element_blank(),
        legend.key.size = unit(2,"line")) + 
  scale_fill_manual(values=c("grey", "grey"), name=NULL) + 
  scale_color_manual(values=c("black", "black")) +
  scale_linetype_manual(values=c("solid", "dashed"), labels = c("AfD vs. Other", "Other vs. AfD")) +
  theme(legend.position = c(0.94, 0.92),
        legend.justification=c(1, 1), 
        legend.background = element_rect(colour="grey80"),
        legend.margin=margin(-3,5,5,5),
        legend.key.size = unit(1.3, "lines")) +
  theme(legend.title=element_blank()) +
  guides(color=FALSE)
plot_party_afd_main_eco_sal_w2

ggsave("./Figures/party_afd_main_distress_w2.pdf", height = 3 , width = 5, scale=1.4, device = cairo_pdf)


#robustness check
plot_party_afd_main_eco_sal_w2_w1 <- ggplot(party_afd_main_eco_sal_w2_w1, aes(x=x, y=predicted, group=group)) +
  theme_minimal() + 
  geom_line(aes(linetype=group, color="black"), size=0.7) +
  geom_ribbon(aes(ymin=party_afd_main_eco_sal_w2_w1$conf.low, ymax=party_afd_main_eco_sal_w2_w1$conf.high), linetype=0, alpha=0.1, show.legend = F) +
  scale_y_continuous(limit = c(-50, 50)) +
  scale_x_continuous(limit = c(1,5), breaks=c(1,2,3,4,5)) +  
  xlab("Economic Distress") + 
  ylab("Sympathy") + 
  labs(color  = "Party", linetype = "Party", fill = "Party", group="Party") +
  theme(legend.position = "bottom", legend.margin=margin(0,0,0,0), legend.box.margin=margin(-5,-5,-5,-5), legend.title=element_blank(),
        legend.key.size = unit(2,"line")) + 
  scale_fill_manual(values=c("grey", "grey"), name=NULL) + 
  scale_color_manual(values=c("black", "black")) +
  scale_linetype_manual(values=c("solid", "dashed"), labels = c("AfD vs. Other", "Other vs. AfD")) +
  theme(legend.position = c(0.94, 0.92),
        legend.justification=c(1, 1), 
        legend.background = element_rect(colour="grey80"),
        legend.margin=margin(-3,5,5,5),
        legend.key.size = unit(1.3, "lines")) +
  theme(legend.title=element_blank()) +
  guides(color=FALSE)
plot_party_afd_main_eco_sal_w2_w1

ggsave("./Figures/party_afd_main_distress_w2_w1.pdf", height = 3 , width = 5, scale=1.4, device = cairo_pdf)








### combine plots
#combine new plots by IV

p_eco_sal <- plot_grid(plot_nat_eco_sal_w1, plot_nat_eco_sal_w2, plot_party_afd_main_eco_sal_w1, plot_party_afd_main_eco_sal_w2,
                       ncol = 2, rel_widths = c(1,1), vjust = 0, labels = c("Wave1", "Wave2"))
p_eco_sal + theme(plot.margin = unit(c(1,0,0,0), "cm"))
ggsave("./Figures/symp_distress.pdf", height = 3 , width = 4.5, scale=2.0, device = cairo_pdf)
ggsave("./Figures/symp_distress.png", height = 3 , width = 4.5, scale=2.0, device = cairo_pdf)



# create Tables

# W1 and W2
texreg(list("Italians" = nat_ita_eco_sal_w1, "Turks" = nat_tur_eco_sal_w1, 
            "Italians" = nat_ita_eco_sal_w2, "Turks" = nat_tur_eco_sal_w2, 
            "AfD vs. Other" = afd_main_eco_sal_w1, "Other vs. AfD" = main_afd_eco_sal_w1,
            "AfD vs. Other" = afd_main_eco_sal_w2, "Other vs. AfD" = main_afd_eco_sal_w2),
       summaries = c("Observations", "R2"), 
       single.row=TRUE, digits = 3, dcolumn = T, caption = "Regression Models of Economic Distress on Sympathy",
       file = "./Tables/symp_distress.tex")

#W1 Effects on W2
texreg(list("Italians" = nat_ita_eco_sal_w2_w1, "Turks" = nat_tur_eco_sal_w2_w1, 
            "AfD vs. Other" = afd_main_eco_sal_w2_w1, "Other vs. AfD" = main_afd_eco_sal_w2_w1),
       summaries = c("Observations", "R2"), 
       single.row=TRUE, digits = 3, dcolumn = T, caption = "Regression Models of Economic Distress on Sympathy",
       file = "./Tables/symp_robustness_check.tex")





#############################################
#############################################
#### Comparisons between 2017 and 2020 ######
#############################################
#############################################


rm(list=ls(all=TRUE))


data <- read_stata("./Data/corona_comp_1720.dta")

de <- data %>% filter(time==1)
d1 <- data %>% filter(time==2)
d2 <- data %>% filter(time==3)


#same nationality
comp_all_de <- lm(feel ~ 1, data = subset(de, group==10))

comp_all <- plot_models(comp_all_de,  
            m.labels = c("DE: 2017"), 
            show.intercept = TRUE, axis.lim=c(-50,50))
comp_all[["data"]]$treatment <- "Same Nationality"
comp_all[["data"]]$n <- c(length(comp_all_de[["model"]][["feel"]]))

#italian
comp_ita_d1 <- lm(feel ~ 1, data = subset(d1, group==9))
comp_ita_d2 <- lm(feel ~ 1, data = subset(d2, group==9))
comp_ita_de <- lm(feel ~ 1, data = subset(de, group==9))
comp_ita <- plot_models(comp_ita_d1, comp_ita_d2, comp_ita_de,  
                        m.labels = c("DE: 2020-W1", "DE: 2020-W2", "DE: 2017"), 
                        show.intercept = TRUE, axis.lim=c(-50,50))
comp_ita[["data"]]$treatment <- "Italian"
comp_ita[["data"]]$n <- c(length(comp_ita_d1[["model"]][["feel"]]),
                          length(comp_ita_d2[["model"]][["feel"]]),
                               length(comp_ita_de[["model"]][["feel"]]))
#turkish
comp_tur_d1 <- lm(feel ~ 1, data = subset(d1, group==8))
comp_tur_d2 <- lm(feel ~ 1, data = subset(d2, group==8))
comp_tur_de <- lm(feel ~ 1, data = subset(de, group==8))
comp_tur <- plot_models(comp_tur_d1, comp_tur_d2, comp_tur_de,  
                        m.labels = c("DE: 2020-W1", "DE: 2020-W2", "DE: 2017"), 
                        show.intercept = TRUE, axis.lim=c(-50,50))
comp_tur[["data"]]$treatment <- "Turkish"
comp_tur[["data"]]$n <- c(length(comp_tur_d1[["model"]][["feel"]]),
                          length(comp_tur_d2[["model"]][["feel"]]),
                               length(comp_tur_de[["model"]][["feel"]]))

#same party
comp_party_equal_d1 <- lm(feel ~ 1, data = subset(d1, ((group==1 & vote==6) | (group==2 & vote==1) | (group==3 & vote==2) | 
                                                         (group==4 & vote==5) | (group==5 & vote==3))))
comp_party_equal_d2 <- lm(feel ~ 1, data = subset(d2, ((group==1 & vote==6) | (group==2 & vote==1) | (group==3 & vote==2) | 
                                                        (group==4 & vote==5) | (group==5 & vote==3))))
comp_party_equal_de <- lm(feel ~ 1, data = subset(de, ((group==1 & vote==6) | (group==2 & vote==1) | (group==3 & vote==2) | 
                                                        (group==4 & vote==5) )))

comp_party_equal <- plot_models(comp_party_equal_d1, comp_party_equal_d2, comp_party_equal_de,  
                        m.labels = c("DE: 2020-W1", "DE: 2020-W2", "DE: 2017", "AT: 2018"), 
                        show.intercept = TRUE, axis.lim=c(-50,50))
comp_party_equal[["data"]]$treatment <- "Same Party"
comp_party_equal[["data"]]$n <- c(length(comp_party_equal_d1[["model"]][["feel"]]),
                                  length(comp_party_equal_d2[["model"]][["feel"]]),
                          length(comp_party_equal_de[["model"]][["feel"]]))

#different party
comp_party_diff_d1 <- lm(feel ~ 1, data = subset(d1, ((group==1 & vote!=6) | (group==2 & vote!=1) | (group==3 & vote!=2) | 
                                                        (group==4 & vote!=5) | (group==5 & vote!=3))))
comp_party_diff_d2 <- lm(feel ~ 1, data = subset(d2, ((group==1 & vote!=6) | (group==2 & vote!=1) | (group==3 & vote!=2) | 
                                                       (group==4 & vote!=5) | (group==5 & vote!=3))))
comp_party_diff_de <- lm(feel ~ 1, data = subset(de, ((group==1 & vote!=6) | (group==2 & vote!=1) | (group==3 & vote!=2) | 
                                                          (group==4 & vote!=5) )))

comp_party_diff <- plot_models(comp_party_diff_d1, comp_party_diff_d2, comp_party_diff_de,  
                                m.labels = c("DE: 2020-W1", "DE: 2020-W2",  "DE: 2017", "AT: 2018"), 
                                show.intercept = TRUE, axis.lim=c(-50,50))
comp_party_diff[["data"]]$treatment <- "Different Party"
comp_party_diff[["data"]]$n <- c(length(comp_party_diff_d1[["model"]][["feel"]]),
                                 length(comp_party_diff_d2[["model"]][["feel"]]),
                          length(comp_party_diff_de[["model"]][["feel"]]))






#############################
## Asymmetric relationship
#############################

#AfD and Mainstream
comp_afd_main_d1 <- lm(feel ~ 1, data = subset(d1, ((group==2 | group==3) & (vote==6)))) 
comp_afd_main_d2 <- lm(feel ~ 1, data = subset(d2, ((group==2 | group==3) & (vote==6)))) 
comp_afd_main_de <- lm(feel ~ 1, data = subset(de, ((group==2 | group==3) & (vote==6)))) 
comp_afd_main <- plot_models(comp_afd_main_d1, comp_afd_main_d2, comp_afd_main_de,  
                             m.labels = c("DE: 2020-W1", "DE: 2020-W2", "DE: 2017"), 
                             show.intercept = TRUE, axis.lim=c(-50,50))
comp_afd_main[["data"]]$treatment <- "AfD vs. Main"
comp_afd_main[["data"]]$n <- c(length(comp_afd_main_d1[["model"]][["feel"]]),
                               length(comp_afd_main_d2[["model"]][["feel"]]),
                               length(comp_afd_main_de[["model"]][["feel"]]))
  
comp_main_afd_d1 <- lm(feel ~ 1, data = subset(d1, (group==1 & (vote==1 | vote==2))))
comp_main_afd_d2 <- lm(feel ~ 1, data = subset(d2, (group==1 & (vote==1 | vote==2))))
comp_main_afd_de <- lm(feel ~ 1, data = subset(de, (group==1 & (vote==1 | vote==2)))) 
comp_main_afd <- plot_models(comp_main_afd_d1, comp_main_afd_d2, comp_main_afd_de,  
                                     m.labels = c("DE: 2020-W1", "DE: 2020-W2", "DE: 2017"), 
                                     show.intercept = TRUE, axis.lim=c(-50,50))
comp_main_afd[["data"]]$treatment <- "Main vs. AfD"
comp_main_afd[["data"]]$n <- c(length(comp_main_afd_d1[["model"]][["feel"]]),
                               length(comp_main_afd_d2[["model"]][["feel"]]),
                               length(comp_main_afd_de[["model"]][["feel"]]))



#AfD and Left-Wing Parties (LP+Greens)
comp_afd_lw_d1 <- lm(feel ~ 1, data = subset(d1, ((group==4 | group==5) & vote==6)))
comp_afd_lw_d2 <- lm(feel ~ 1, data = subset(d2, ((group==4 | group==5) & vote==6)))
comp_afd_lw_de <- lm(feel ~ 1, data = subset(de, ((group==4 | group==5) & vote==6))) 
comp_afd_lw <- plot_models(comp_afd_lw_d1, comp_afd_lw_d2, comp_afd_lw_de,  
                             m.labels = c("DE: 2020-W1", "DE: 2020-W2", "DE: 2017"), 
                             show.intercept = TRUE, axis.lim=c(-50,50))
comp_afd_lw[["data"]]$treatment <- "AfD vs. Left-Wing"
comp_afd_lw[["data"]]$n <- c(length(comp_afd_lw_d1[["model"]][["feel"]]),
                               length(comp_afd_lw_d2[["model"]][["feel"]]),
                               length(comp_afd_lw_de[["model"]][["feel"]]))

comp_lw_afd_d1 <- lm(feel ~ 1, data = subset(d1, (group==1 & (vote==3 | vote==5))))
comp_lw_afd_d2 <- lm(feel ~ 1, data = subset(d2, (group==1 & (vote==3 | vote==5))))
comp_lw_afd_de <- lm(feel ~ 1, data = subset(de, (group==1 & (vote==3 | vote==5)))) 
comp_lw_afd <- plot_models(comp_lw_afd_d1, comp_lw_afd_d2,  comp_lw_afd_de,  
                             m.labels = c("DE: 2020-W1", "DE: 2020-W2", "DE: 2017"), 
                             show.intercept = TRUE, axis.lim=c(-50,50))
comp_lw_afd[["data"]]$treatment <- "Left-Wing vs. AfD"
comp_lw_afd[["data"]]$n <- c(length(comp_lw_afd_d1[["model"]][["feel"]]),
                               length(comp_lw_afd_d2[["model"]][["feel"]]),
                               length(comp_lw_afd_de[["model"]][["feel"]]))


### aggregate results
comp_results_asym_v2 <- data.table::rbindlist(list(comp_all[["data"]], comp_ita[["data"]], comp_tur[["data"]], comp_party_equal[["data"]], comp_party_diff[["data"]],
                              comp_afd_main[["data"]], comp_main_afd[["data"]], comp_afd_lw[["data"]], comp_lw_afd[["data"]]), fill=T)

comp_results_asym_v2_d1 <- comp_results_asym_v2 %>% filter(group == "DE: 2020-W1")
comp_results_asym_v2_d2 <- comp_results_asym_v2 %>% filter(group == "DE: 2020-W2")
comp_results_asym_v2_de <- comp_results_asym_v2 %>% filter(group == "DE: 2017")

comp_results_asym_v2_d1$term <- comp_results_asym_v2_d1$treatment
comp_results_asym_v2_d2$term <- comp_results_asym_v2_d2$treatment
comp_results_asym_v2_de$term <- comp_results_asym_v2_de$treatment

comp_results_asym_v2_d1$model <- "DE: 2020-W1"
comp_results_asym_v2_d2$model <- "DE: 2020-W2"
comp_results_asym_v2_de$model <- "DE: 2017"



# change labels
comp_asym_models_v2_de <- data.table::rbindlist(list(comp_results_asym_v2_d1, comp_results_asym_v2_d2, comp_results_asym_v2_de), fill=T)
comp_asym_models_v2_de$term <- recode(comp_asym_models_v2_de$term, "Same Nationality" = "Nationality: Same", 
                                   "Italian" = "Italian", "Turkish" = "Turkish",
                                   "Same Party" = "Party: Same", "Different Party" = "Different", 
                                   "AfD vs. Main" = "AfD towards Major Parties", "Main vs. AfD" = "Major Parties towards AfD",
                                   "AfD vs. Left-Wing" = "AfD towards Left-Wing", "Left-Wing vs. AfD" = "Left-Wing towards AfD")

comp_asym_models_v2_de <- comp_asym_models_v2_de %>% add_row(estimate = -99, term = "Nationality: Same", model = "DE: 2020-W1", group = "DE: 2020-W1")
comp_asym_models_v2_de <- comp_asym_models_v2_de %>% add_row(estimate = -99, term = "Nationality: Same", model = "DE: 2020-W2", group = "DE: 2020-W2")

comp_asym_models_v3_de <- rbind(comp_asym_models_v2_de[17:25,],comp_asym_models_v2_de[26,], 
                                comp_asym_models_v2_de[1:8,],comp_asym_models_v2_de[27,], comp_asym_models_v2_de[9:16,])

comp_asym_models_v3_de$model <- factor(comp_asym_models_v3_de$group)


#create plot
dwplot(comp_asym_models_v3_de, dot_args = list(aes(shape = model), color="black", fill="white", size=2),
       whisker_args = list(aes(linetype = model), color="black"), dodge_size = 0.6, vline = geom_vline(xintercept=0, color = "grey60")) +
  xlab("Feeling Thermometer") + 
  scale_x_continuous(limits = c(-50,50)) +
  scale_linetype_manual(values=c("solid", "dashed", "dotted"), 
                        guide = guide_legend(reverse = TRUE)) +
  theme_minimal() + 
  theme(legend.position = c(0.945, 0.2),
        legend.justification=c(1, 1), 
        legend.background = element_rect(colour="grey80"),
        legend.margin=margin(-3,5,5,5),
        legend.key.size = unit(1, 'lines')) +
  theme(legend.title=element_blank()) 

ggsave("./Figures/time_comparison_2017-20.pdf", height = 2 , width = 3, scale=2.2, device = cairo_pdf)
ggsave("./Figures/time_comparison_2017-20.png", height = 2 , width = 3, scale=2.2, device = cairo_pdf)




######################
### create Tables ####
######################


df1 <- d1 %>% select(Vignette, edu2, sex, age, income, cducsu, spd, fdp, green, left, afd, other, dnv, feel)
df2 <- d2 %>% select(Vignette, edu2, sex, age, income, cducsu, spd, fdp, green, left, afd, other, dnv, feel)
df17<- de %>% select(Vignette, edu2, sex, age, income, cducsu, spd, fdp, green, left, afd, other, dnv, feel)


#############################
# time comparison models ####
#############################

export_df <- comp_asym_models_v3_de %>% mutate_at(vars(estimate, std.error, conf.low, conf.high), funs(round(., 3)))

export_d2 <- export_df %>% filter(group=="DE: 2020-W2") %>% select(treatment, estimate, std.error, conf.low, conf.high, n)
export_d1 <- export_df %>% filter(group=="DE: 2020-W1") %>% select(treatment, estimate, std.error, conf.low, conf.high, n)
export_de <- export_df %>% filter(group=="DE: 2017") %>% select(treatment, estimate, std.error, conf.low, conf.high, n)




latex(export_d2, file="./Tables/comp_models_2020-W2.tex", 
      rowname = NULL, dcolumn=T, caption = "Regression Models 2020: Wave 2", label = "tab:comp_models_2020_w2", where="!htbp")
latex(export_d1, file="./Tables/comp_models_2020-W1.tex", 
      rowname = NULL, dcolumn=T, caption = "Regression Models 2020: Wave 1", label = "tab:comp_models_2020_w1", where="!htbp")
latex(export_de, file="./Tables/comp_models_2017.tex", 
      rowname = NULL, dcolumn=T, caption = "Regression Models 2017", label = "tab:comp_models_2017", where="!htbp")


########################
#summary statistics ####
########################

df2 %>%  select(sex, age, edu2, income, feel) %>% 
  stargazer(out="./Tables/summaries_2020-W2.tex",
            align = T)
df1 %>% select(sex, age, edu2, income, feel) %>% 
  stargazer(out="./Tables/summaries_2020-W1.tex",
            align = T)
de %>% select(sex, age, edu2, income, feel) %>% 
  stargazer(out="./Tables/summaries_2017.tex",
            align = T)


##########################
#treatment assignment ####
##########################

#2020: W2
treat_socdem_d2 <- df2 %>%                        # Specify data frame
  group_by(Vignette) %>%                         # Specify group indicator
  summarise_at(vars(edu2, sex, age, cducsu, spd, fdp, green, left, afd, other, dnv), c(mean), na.rm = TRUE)   

colnames(treat_socdem_d2) <- c("Treatment", "Education", "Sex", "Age", "CDUCSU", "SPD", "FDP", "Greens", "LeftParty", "AfD", "Other", "donotvote")

treat_socdem_d2$Treatment <- c("City/AfD", "City/CDU", "City/SPD", "City/Left Party", "City/Greens", "City/EU", "City/never foreign country", 
                              "Rural/AfD", "Rural/CDU", "Rural/SPD", "Rural/Left Party", "Rural/Greens", "Rural/EU", "Rural/never foreign country",
                              "Turkey", "Italy")

treat_socdem_d2 <- treat_socdem_d2 %>% mutate_at(vars(Sex, Age, Education, CDUCSU, SPD, FDP, Greens, LeftParty, AfD, Other, donotvote), funs(round(., 3)))

latex(treat_socdem_d2, file="./Tables/treat_socdem_2020-W2.tex", 
      rowname = NULL, dcolumn=T, caption = "Treatment Assignment by Socio-Demographics 2020 (Wave 2)", label = "tab:treat_socdem_2020_W2", where="!htbp")

#2020: W1
treat_socdem_d1 <- df1 %>%                        # Specify data frame
  group_by(Vignette) %>%                         # Specify group indicator
  summarise_at(vars(edu2, sex, age, cducsu, spd, fdp, green, left, afd, other, dnv), c(mean), na.rm = TRUE)   

colnames(treat_socdem_d1) <- c("Treatment", "Education", "Sex", "Age", "CDUCSU", "SPD", "FDP", "Greens", "LeftParty", "AfD", "Other", "donotvote")

treat_socdem_d1$Treatment <- c("City/AfD", "City/CDU", "City/SPD", "City/Left Party", "City/Greens", "City/EU", "City/never foreign country", 
                               "Rural/AfD", "Rural/CDU", "Rural/SPD", "Rural/Left Party", "Rural/Greens", "Rural/EU", "Rural/never foreign country",
                               "Turkey", "Italy")

treat_socdem_d1 <- treat_socdem_d1 %>% mutate_at(vars(Sex, Age, Education, CDUCSU, SPD, FDP, Greens, LeftParty, AfD, Other, donotvote), funs(round(., 3)))

latex(treat_socdem_d1, file="./Tables/treat_socdem_2020-W1.tex", 
      rowname = NULL, dcolumn=T, caption = "Treatment Assignment by Socio-Demographics 2020 (Wave 1)", label = "tab:treat_socdem_2020_W1", where="!htbp")

#2017
treat_socdem_de <- df17 %>%                        # Specify data frame
  group_by(Vignette) %>%                         # Specify group indicator
  summarise_at(vars(edu2, sex, age, cducsu, spd, fdp, green, left, afd, other, dnv), c(mean), na.rm = TRUE)   

colnames(treat_socdem_de) <- c("Treatment", "Education", "Sex", "Age", "CDUCSU", "SPD", "FDP", "Greens", "LeftParty", "AfD", "Other", "donotvote")

treat_socdem_de <- treat_socdem_de %>% mutate_at(vars(Sex, Age, Education, CDUCSU, SPD, FDP, Greens, LeftParty, AfD, Other, donotvote), funs(round(., 3)))

latex(treat_socdem_de, file="./Tables/treat_socdem_de.tex", 
      rowname = NULL, dcolumn=T, caption = "Treatment Assignment by Socio-Demographics 2017", label = "tab:treat_socdem_2017", where="!htbp")





